


first.difference.table <- function(model.list, low.val, high.val){ 
     
     max = 0
     
     for(i in 1: length(model.list))
          if(ncol(model.list[[i]]$model)-1 > max )
               max = ncol(model.list[[i]]$model)-2 

          first.dif  = matrix(NA, nrow = max, ncol = length(model.list))
          
          for(i in 1: length(model.list)){  
               
               some.data = model.list[[i]]$model
               
               for(j in 2: ncol(some.data)-1){ 
                    
                    data.low = some.data
                    data.high = some.data
                    
                    data.low[,j] = quantile(some.data[,j] , probs = low.val, na.rm=T, names= F)
                    data.high[,j]  = quantile(some.data[,j] , probs = high.val, na.rm=T, names= F)
                    
                    pred.low = predict(model.list[[i]], newdata = data.low, type = "response", se.fit = F)
                    pred.high = predict(model.list[[i]], newdata = data.high, type = "response", se.fit = F)
                    
                    first.dif[j-1,i] = (mean(pred.high) - mean(pred.low) ) * 100
               }
          }
          
          first.dif = round(first.dif, 1)
          first.dif[is.na(first.dif)] = "-"
          
          star.table  = matrix(1, nrow = nrow(first.dif), ncol = length(model.list))
          for(i in 1: length(model.list))
               star.table[1:length(coef(summary(model.list[[i]]))[-1,4]),i] = coef(summary(model.list[[i]]))[-1,4]
          
          for(i in 1:nrow(first.dif))
               for(j in 1:ncol(first.dif))
                    if(star.table[i,j] <=.05)
                         first.dif[i,j] = paste0("startbold", first.dif[i,j], "endbold")
          
          return(first.dif)
}




format.first.diff <- function(table, title, models){
     
     lastline = paste0("\\multicolumn{", models+1, "}{l}{Note: Values in bold are significant at $\\alpha \\leq .05$}\\\\" )
     table[1:2] = ""
     table[3] = "\\begin{table}[ht]"
     table[5] = paste0(paste0("\\caption{",title,"}    \\label{}   ", table[5]))
     table[5] = sub(pattern = "\\\\begin\\{tabular\\}\\{", x=table[5], replacement = "\\\\begin\\{tabular\\}\\{@\\{\\\\extracolsep\\{5pt\\}\\}")
     table[6] = "\\toprule "
     table[8] = "\\midrule"
     table[length(table)-2] = paste0(" \\hline \\hline \\\\[-1.8ex]", lastline)
     table = gsub(pattern = "\\$\\\\backslash\\$textbf\\\\", x = table, replacement = "\\\\textbf")
     table = gsub(pattern = "\\\\\\}", x = table, replacement = "\\}")
     table = gsub(pattern = "startbold", x = table, replacement = "\\\\textbf\\{")
     table = gsub(pattern = "endbold", x = table, replacement = "\\}")
     cat(table)
}


change.names <- function(some.model, coef.names){
     names(some.model$coefficients) = coef.names
     return(some.model)
}



pseudo.R.quasi <- function(object){
     #update the function; call as a nonquasi function
     new.model = update(object, family = binomial(link="logit"))
     
     # calculate the dispersion parameter from a non quasi model
     d = with(new.model, sum((weights * residuals^2)[weights > 0])/df.residual)
     
     #calcualte the null 
     new.null = glm(new.model$model[,1] ~ 1, family = binomial(link="logit"), weights = new.model$prior.weights)    
     
     #calculate the quasi AIC
     quasi.AIC.model = QAIC(new.model, chat = d)
     quasi.AIC.null = QAIC(new.null, chat = d)
     
     LL.model = (quasi.AIC.model - 2*length(new.model$coefficients))/-2
     LL.null = (quasi.AIC.null - 2*length(new.null$coefficients))/-2
     
     return(1 - (LL.model/LL.null))
}


stars <- function(x){
     if (x > .1)
          return('')
     else if (x > .05)
          return('*')
     else if (x >.01)
          return('**')
     else 
          return('***')
     
}
